library(tidyverse)
library(tidymodels)
library(modelr)
library(ISLR)
library(GGally)
library(pROC)
library(cowplot)
library(OneR)
library(rlang)
library(caret)
set.seed(1992)
La regresión logística es útil para problemas de predicción de clases. En este caso vamos a tratar de resolver el problema de predecir si una persona va defaultear su deuda de tarjeta de crédito en base a ciertos predictores.
Este conjunto de datos proviene de la librería ISLR (Introduction to Statistical Learning Using R) de James, Witten, Hastie y Tibshirani.
# cargamos los datos y observamos su estructura
default <- Default
glimpse(default)
Rows: 10,000
Columns: 4
$ default [3m[38;5;246m<fct>[39m[23m No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, …
$ student [3m[38;5;246m<fct>[39m[23m No, Yes, No, No, No, Yes, No, Yes, No, No, Yes, Yes, No, No, No, No, No, Yes, No, No, No, No, No, No, No, No, No…
$ balance [3m[38;5;246m<dbl>[39m[23m 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 919.5885, 825.5133, 808.6675, 1161.0579, 0.0000, 0.0000, 1220…
$ income [3m[38;5;246m<dbl>[39m[23m 44361.625, 12106.135, 31767.139, 35704.494, 38463.496, 7491.559, 24905.227, 17600.451, 37468.529, 29275.268, 218…
Tiene 4 variables:
default: La clase que queremos predecir.
student: Binaria que indica si la persona es estudiante.
balance: Balance promedio que le queda a la persona luego de sus pagos mensuales.
income: Ingreso de la persona.
Analicemos la distribución de la clase.
default %>%
group_by(default) %>%
summarise(numero_casos=n())
`summarise()` ungrouping output (override with `.groups` argument)
Vemos que estamos trabajando con un problema de clasificación con un claro desbalance de clase.
Realizamos un gráfico exploratorio completo para ver el comportamiento y las relaciones entre las variables. El color rojo designa a quienes no defaultean y el azul a los que sí.
# graficamos con ggpairs coloreando por variable a predecir
g <- ggpairs(default, 1:4, title = "Correlograma de variables",
mapping = aes(colour= default)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
theme_bw()
# hacemos un loop para cambiar los colores del gráfico
for(i in 1:g$nrow) {
for(j in 1:g$ncol){
g[i,j] <- g[i,j] +
scale_fill_brewer(palette="Set1") +
scale_color_brewer(palette="Set1")
}
}
g
¿Qué pueden decir de la relación entre balance y default?
¿Y entre income y balance?
¿Cuáles parecen ser buenas variables para discriminar entre quienes defaultean y quienes no?
Para modelizar va a ser necesario tener la variable default como numérica. Definimos la variable como {0,1} para los valores {“No”,“Yes”}
default <- default %>%
mutate(default = case_when(default == "No" ~ 0,
default == "Yes" ~ 1))
default
Queremos estimar \(P(Default=Yes|X)=P(X)\) para cada individuo y a partir de ello poder definir un punto de corte para predecir quiénes son los que van a entrar en default.
En este caso estamos modelando la probabilidad de la siguiente manera:
\(P(X)= \beta_0 + \sum\limits_{j=1}^p \beta_j X\)
Veamos que tan bueno es el modelo lineal para esto, usando balance como predictor.
test_mco <- default %>%
lm(formula = default ~ balance, data = .)
tdy = test_mco %>% tidy()
tdy
test_mco %>% glance()
Ambos estimadores son significativos y el test de significatividad global del modelo también es significativo.
Veamos un gráfico de nuestro modelo.
Parece tener bastantes problemas para estimar la probabilidad de default de los individuos. Por ejemplo, vemos que hay varios individuos a los cuales les asigna una probabilidad negativa.
Para evitar estos problemas, usamos la función logística.
\(P(Y=1|X)= \frac{e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}{1+e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}\)
El lado derecho se llama expit
Esta función acota el resultado entre 0 y 1, lo cual es mucho más adecuado para modelar una probabilidad.
Luego de hacer algunas operaciones, podemos llegar a la expresión:
\(\log {\frac{P(x)}{1-P(x)}}= \beta_0 + \sum\limits_{j=1}^p \beta_j X\)
El lado izquierdo es el logaritmo de los odds y se llama logit.
Realizamos una partición una partición entre dataset de entrenamiento (70%) y testeo (30%) usando la función initial_split del paquete rsample de tidymodels.
# fijamos semilla
set.seed(44)
# Partición Train y Test, indicando proporción
train_test <- initial_split(default, prop = 0.7)
# armamos dataframe de testeo y entrenamiento
default <- training(train_test)
test <- testing(train_test)
# vemos el contenido
default %>%
dim_desc() # 7000 filas
[1] "[7,000 x 4]"
test %>%
dim_desc() # 3000 filas
[1] "[3,000 x 4]"
Para aplicar la regresión logística primero usamos la función formulas del paquete modelr para crear un objeto que contiene todas las fórmulas que vamos a utilizar. En .response especificamos la variable respuesta de nuestras fórmulas y luego nombramos las fórmulas que queramos armar.
Así, armaremos distintos modelos combinando distintas variables.
logit_formulas <- formulas(.response = ~ default, # único lado derecho de las formulas.
bal = ~ balance,
stud = ~ student,
inc = ~ income,
bal_stud = ~ balance + student,
bal_inc = ~ balance + income,
stud_inc = ~ student + income,
full = ~ balance + income + student
)
logit_formulas # observamos el objeto formulas
$bal
default ~ balance
$stud
default ~ student
$inc
default ~ income
$bal_stud
default ~ balance + student
$bal_inc
default ~ balance + income
$stud_inc
default ~ student + income
$full
default ~ balance + income + student
Procedemos a crear los modelos a partir de estas fórmulas.
models <- data_frame(logit_formulas) %>% # dataframe a partir del objeto formulas
mutate(models = names(logit_formulas), # columna con los nombres de las formulas
expression = paste(logit_formulas), # columna con las expresiones de las formulas
mod = map(logit_formulas, ~glm(., family = 'binomial', data = default))) # Que estamos haciendo acá? Que vamos a encontrar en la columna?
models
La funcíón glm() nos permite crear un modelo lineal generalizado (Generalized Linear Model). Al igual que la función lm() toma como argumentos una formula y los datos pero también se debe especificar el argumento family: indicamos la distribución del error y la función link que vamos a utilizar en el modelo.
Algunas familias son:
Binomial: link=logit
Poisson: link=log
Gaussiana: link=identidad
Como estamos trabajando con un fenómeno que suponemos tiene una distribución binomial, así lo especificamos en el parámetro family.
Probamos los primeros tres modelos, aquellos que tienen un único predictor. Usamos la función tidy para obtener los parámetros estimados para estos tres modelos.
models %>%
filter(models %in% c('bal','stud','inc')) %>%
mutate(tidy = map(mod, tidy)) # Qué realizamos en este paso? Que va a tener esta columna?
Para acceder a los elementos de la nueva columna tidy debemos desanidarla.
models %>%
filter(models %in% c('bal','stud','inc')) %>%
mutate(tidy = map(mod, tidy)) %>% # Qué realizamos en este paso? Que va a tener esta columna?
unnest(tidy) %>%
mutate(estimate=round(estimate,5), # redondeamos valores para facilitar lectura
p.value=round(p.value,4))
Observamos que todos los modelos tienen coeficientes significativos excepto el coeficiente asociado al Ingreso (income), el cual tiene un p-valor de 0.17.
Recordando la ecuación para modelar la probabilidad:
\(P(Y=1|X)= \frac{e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}{1+e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}\)
Se observa que ahora las variables ya no tienen una relación lineal con la probabilidad. En este modelo un coeficiente positivo indica que frente a aumentos de dicha variable la probabilidad aumenta, mientras que un coeficiente negativo nos indica lo contrario. Para nuestro caso:
El coeficiente de 0.0056 de la variable balance indica que cuando el balance aumenta la probabilidad de default aumenta.
El coeficiente de 0.40886 de la variable studentYes indica que cuando una persona es estudiante su probabilidad de default aumenta respecto a una persona que no lo es.
El coeficiente de -0.00001 de la variable income indica que cuando el ingreso aumenta la probabilidad de default disminuye.
Recomendamos leer el capítulo de regresión logística de Interpretable Machine Learning: A Guide for Making Black Box Models Explainable de Molnar Cristoph para una discusión más profunda de la interpretación de un modelo de regresión logística.
Ahora probamos con un modelo que utiliza las tres predictoras.
models %>%
filter(models == "full") %>%
mutate(tidy = map(mod,tidy)) %>%
unnest(tidy) %>%
mutate(estimate=round(estimate,5),
p.value=round(p.value,4))